#!/usr/bin/perl
# Sep 29, 2009 Overlap analysis

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
use Math::Complex;

my $search_width = 300;


if((@ARGV != 4 and @ARGV != 5) or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [bed file #1] -b [bed file #2] -q (silent mode)\n";
}

my %opt;
getopts("a:b:q", \%opt);

my %chr_conv = (
	'I'	=>	'chr1',
	'II'	=>	'chr2',
	'III'	=>	'chr3'
);


#---------------------------------------
# region array making
#---------------------------------------
#chr1    0       200     200     1239.11 4.358   1.289999e-114   2.02396394827586e-114
#chr1    1910    2090    1910    1122.16 4.318   1.133638e-89    1.55519182914573e-89

my %all;
my %regionA = &readBed($opt{a});
my %regionB = &readBed($opt{b});
sub readBed{
	my ($file) = @_;
	my %hash;
	my $fh = IO::File->new($file) or die "cannot open $file: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/^track/){
			next;
		}
		my ($chr, $start, $end, @other) = split /\t/;
		if(exists $chr_conv{$chr}){
			$chr = $chr_conv{$chr};
		}

		$all{$chr}{$start - 150}++;
		$all{$chr}{$end + 149}++;
		$hash{$chr}{$start - 150} = 0;	#start
		$hash{$chr}{$end + 149} = 1;		#end
	}
	$fh->close();
	return %hash;
}


#---------------------------------------
# merge peaks
#---------------------------------------
my $common = 0;
my $uniqueA = 0;
my $uniqueB = 0;
my $flag_inA = 0;
my $flag_inB = 0;
my $flag_common = 0;
my $start = 0;
my $end = 0;
my %occupy;
my $SUM_peaks = 0;
my $NUM_peaks = 0;
foreach my $chr(qw(chr1 chr2 chr3)){
	foreach my $key (sort {$a <=> $b} keys %{$all{$chr}}){
		if(exists $regionA{$chr}{$key} and $regionA{$chr}{$key} == 0){
			$flag_inA = 1;
			if($flag_inB == 1){
				$flag_common = 1;
			}else{
				$start = $key;
			}
		}
		if(exists $regionB{$chr}{$key} and $regionB{$chr}{$key} == 0){
			$flag_inB = 1;
			if($flag_inA == 1){
				$flag_common = 1;
			}else{
				$start = $key;
			}
		}
		if(exists $regionA{$chr}{$key} and $regionA{$chr}{$key} == 1){
			$flag_inA = 0;
			if($flag_common == 0){
				$uniqueA++;
				$end = $key;
			}elsif($flag_inB == 1){
				# 共通でもう１方はまだ続いている
			}elsif($flag_inB == 0){
				# 共通でもう１方はすでに終わった
				$common++;
				$flag_common = 0;
				$end = $key;
			}
		}
		if(exists $regionB{$chr}{$key} and $regionB{$chr}{$key} == 1){
			$flag_inB = 0;
			if($flag_common == 0){
				$uniqueB++;
				$end = $key;
			}elsif($flag_inA == 1){
				# 共通でもう１方はまだ続いている
			}elsif($flag_inA == 0){
				# 共通でもう１方はすでに終わった
				$common++;
				$flag_common = 0;
				$end = $key;
			}
		}
		if($end != 0){
			for(my $i = $start + 150; $i <= $end - 149; $i++){
				$occupy{$chr}{$i} = 1;
			}
			$SUM_peaks += $end - $start - 300;
			$NUM_peaks++;
			$end = 0;
		}
	}
}

#---------------------------------------
# calculation
#---------------------------------------
my $genome_length = 12781358;
my $non_peaks_length = $genome_length - $SUM_peaks;
my $average_peak_length = $SUM_peaks / $NUM_peaks;
my $estimated_non_peak_number = int($non_peaks_length / $average_peak_length);

#---------------------------------------
# output result
#---------------------------------------
my $R_string = 'y$p.value';
my $Rcommand_p = <<"END";
R --vanilla --slave <<'EOF'
 Convictions <- matrix(c($common,$uniqueA,$uniqueB,$estimated_non_peak_number), nrow=2, dimnames=list(c("common", "uniqueA"), c("uniqueB","other")))
 y <-  fisher.test(Convictions, alternative="greater")
 cat($R_string)
EOF
END

my $pvalue = `$Rcommand_p`;

if($opt{q}){
	print "$uniqueA\t$common\t$uniqueB\t$pvalue\n";
}else{
	print "common: $common\n";
	print "uniqueA: $uniqueA\n";
	print "uniqueB: $uniqueB\n";
	print "other: $estimated_non_peak_number\n";
	printf "average peak length: %.2f\n", $average_peak_length;
	print "p-value = $pvalue\n";
}


